The Excursion Set Mass Function of Superclusters 

Seunghwan Lim and Jounghun Lee 
Department of Physics and Astronomy, Seoul National University, Seoul 151-747, Korea 
shlim@astro . snu . ac . kr , j ounghun@astro . snu . ac . kr 

ABSTRACT 

The mass function of superclusters is derived fully analytically with the help 
of the extended excursion set theory and shown to be in excellent agreement with 
the numerical results from various publicly available N-body simulation database. 
We introduce a new multi-dimensional barrier model in which the formation of su- 
perclusters occurs when the initial shear eigenvalues that perform non-Markovian 
random walks enter a clustering zone surrounded by one reflecting and two ab- 
sorbing barriers. The multi-dimensional barrier heights are determined from the 
first order Lagrangian perturbation theory and found to be independent of red- 
shift and background cosmology. With the help of our analytic model for the 
supercluster mass function, the relative abundance of the rich superclusters is 
analytically evaluated at a given epoch and found to be sensitive to the growth 
rate of the cosmic web. Our result implies that the relative abundance of the rich 
superclusters at a given epoch may be useful as a cosmological test of gravity. 

Subject headings: cosmology:theory — large scale structure of Universe 



INTRODUCTION 



The gravitational aggregates of (a few to hundreds of) galaxy clusters are often called the 
superclusters which are marginally bound systems, sources of soft X-ray and the Sunyaev 



Zel'dovich effect (e.g., iMyers et al.l l2004t IZappacosta et al 



2005; 



Sadeh fc Rephaelil 12005 



Einasto et al. 



2011 



and references 



see), and have conspicuously filamentary shapes (e.g., 
therein). Although quite rare in the local Universe, the superclusters are believed to be 
quite common phenomena on the scales larger than 100 Mpc. The nearby Virgo cluster as 
well as the Local Group where our Milky Way r esides also b elongs to the Local Superclusters 
that contains more than 100 member clusters (TuUyj|l982|, and references there in). 



In a ACDM (A+cold dark matter) universe, the formation of the superclusters at the 
present epoch represents the grand finale of the hierarchical merging events. No bound 
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objects could form on mass scales larger than that of the present rich superclusters in the 
future due to the anti-gravitational effect of the cosmological constant (A). That is, the 
rich superclusters observed at present epoch will end up as isolate d massive clusters in the 
future when the A becomes progressively more dominant (e.g., see 
Busha et al]l2005[ krava-Melo et al.ll2009h . 



Nagamine fc Loebl 12003 



If the dark energy were not A or if the large-scale gravity deviates from the general 
relativity (GR), the superclusters could meet a different fate. For example, in QCDM 
(Quintessence+CDM) models, more massive objects than the rich sup erclusters would form 



through the large scale clustering of the Quintessence scalar field (see ICaldwell et al.lll998 



and references therein). In some modified gravity scenarios where the Universe has no anti- 
gravitational dark energy (jMilgroml Il983l ). nothing would prevent the superclusters from 
assembling into larger scale objects. Henceforth, the number counts of the rich superclusters 
might be a powerful indicator of any modification to gravity and the true nature of dark 
energy. 

There are two advantages that the number counts of the rich superclusters have as a 
cosmological probe over that of the clusters. First of all, the rich superclusters are larger and 
rarer on average than the clusters, and thus their abundance should be more sensitive to the 
growth rate and clustering of the large-scale structure. Second, since the superclusters are 
still in the quasi-linear regime, their mass function may be fully analytically tractable with 
the Lagrangian perturbation theory alone. It goes without saying that a robust analytic 
model derived from the first principles is much more desirable than the empirical fitting 
formula. 



It was lOguri et al.l ( 120041 ) who for the first time attempted to analytically derive the 
supercluster mass function (i.e., the number density of the superclusters per unit volume as 
a function of mass), in th e framework of the standard Press-Schechter mass function theory 
(IPress fc Schechterl Il974j ). Pioneering as it was, their work was based on a few unjustified 
assumptions, which led their model to fa il in matching the N-body results (M. Oguri in 
private communication). lYan fc Fan! (120111 ) determined numerically the mass function of the 
supercluster-like filaments using the data from large N-body simulations and showed that the 
standard excursion set theory is incapable of reproducing the numerical results. Although 
they claimed that that when the peak-exclusion effect is incorporated, the excursion set mass 
function of the superclusters shows a qualitative agreement with the numerical results, their 
model still fails in making robust quantitative success. The main difficulty in modeling the 
supercluster mass function lies in the fact that the superclusters are not relaxed systems 
but in the middle of triaxial collapsing process. Instead of just assuming naively that the 
superclusters correspond to those less overdense sites in the linear density field than the 
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clusters, the dynamical state of the on-going triaxial collapse process has to be taken into 
account for the evaluation of the supercluster mass function. 

Here our aim is to develop a robust analytical model for the supercluster mass function 
that would match the results from the largest-scale N-body simulations. To accomplish this 
aim, the modification of the excursion set mass function formalism is made in the direction 
of incorporating the triaxial clustering conditions for the formation of the superclusters with 
the help of the first order Lagrangian theory. This Paper is organized into four sections 
whose contents are summarized as follows: In §5] we present an analytic model for the 
supercluster mass function through modifying the excursion set formalism. In £J3] we test 
the analytic model against the numerical results obtained utilizing the publicly available 
data from various N-body simulations. In §H we explore the possibility of using the relative 
abundance of the rich superclusters as a cosmological test of gravity. In §5] the prospects of 
our work are discussed and the bottom line is stated. 



AN ANALYTIC MODEL 



2.1. Key Concepts of the Excursion Set Formalism 



The standard excursion set mass functio n theory has two key conc e pts: the random- 



walks and the absorbing barrier ([Adlerl Il98ll : iPeacock fc Heavens! Il990t iBond et al.l 11991 



Jedamzikill995l : lRednerll200ll : IZentnerll2007l ) . The growth of the linearly extrapolated density 
contrast, 5(x) = Ap(x)/p, in a given Lagrangian region, x, is described as diffusive random 
walking process as a function of the time-like variable, o~(M), that monotonically increases 
as the random walking proceeds. Here, p(x) is a smoothed density field, p is a mean mass 
density of the Universe, and <r(M) is the rms fluctuation of the linear density, defined as 

1 



a\M) 



2tt 2 



oik k 2 P(k)W 2 (k;M), 



1 



where P{k) is the power spectrum of <5(x) that depends on the key cosmological parameters, 
and W(k; M) is the filter (in Fourier space) with which <5(x) is smoothed on the mass scale 
of M. Throughout this Paper, we choose the t op-hat spherical w indow function for W(k; M) 
and utilize the publicly available CAMB code (ILewis et al.ll2000l ) to compute P(k) and a(M). 



When a random walk hits a prescribed absorbing barrier with height of 5 C at cr(M), 
it corresponds to the occurrence of the gravitational collapse to form a bound object of 
mass M. By saying an absorbing barrier, it means that the barrier represents the lower 
limit of 5: once the value of 5 exceeds the lower limit, 5 > 5 C , the barrier absorbs the 
density growth, bringing its random walking to a halt. The barrier height (i.e., the density 
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threshold), 8 C , is determined by the underlying dynamics of the gravi tational collapse. It 
has a constant value of S r = 1.68 6 in the spherical collapse model (IGunn &: Gottl Il972 



Peebles! 11980 
models (e.g., 



Bardeen et al.l 119861). while it depends on a(M) in the ell i psoidal collapse 



Shethfc Tormenlll999l: fchiueh & LeelhoOll : Isheth et al.lboQll : iLin etaD 12002 



Robertson etal1l2009l ; iLudlow fc Porciani|[20fl] ). 



In this excursion set picture the mass function of bound objects, dN/dM, can be ob- 
tained just by counting the number density of those random walks that hit the absorbing 
barrier in a differential interval of [cr, a + da] . Mathematically, this excursion set mass func- 
tion of bound objects, dN/dM, is a solution to the following integral Equation (IJedamzik 
1995h 

r°° M' dN 

F(5 C ;M)= dM'—^-P(M,M'). (2) 
J M p dM' 

Here F(5 C ; M) is the volume fraction oc cupied by those initial regions with 5 > 5 C on the 
mass scale of M (IPress &: Schechterlll974l ). (M' / p)dN/dM' is the differential volume fraction 
occupied by those regions with 5' = 5 C on some larger mass scale of M' > M, and P{M, M') 
is the conditional probabi lity of having S > 8 r . on the mass scale of M provided that 5' = 5 C 
on the larger scale of M' (jBowerl Il99ll ; lLacey fc Cold I19931 ) : 



F{5 C ;M) 
P(M, M') 



dSp[S;a(M)l 

d5p[5,a(M)\5' = 8 c ,a(M% 



(3) 
(4) 



It is worth recalling that Equation ()2]) formulated first by iJedamzikl ( 119951 ) is com- 
pletely equivalent to the diffusion equation for the mass function of the bound objects 
( iPeacock fc Heavens! Il990l ; iBond et al.l I19911 ) which explicitly shows that the growth of 5 
is a diffusive random walking process. Equation is a reformulation of the diffusion equa- 
tion, depicting implicitly the random walking process of 5. For mathematical convenience, 
we prefer using Jedamzik's formalism, i.e., the implicit version of the diffusion equation, for 
the evaluation of dN/dM. 

The excursion set mass function derived as a solution to Equation ([2]) automatically 
satisfies the normalization condition of J °° dF(8 c , M) = 1, properly taking into account 
the occurrences of the clouds- in- clouds (underdense regions embedded in overdense regions), 
which in fact is reflected by the conditional probability, P(M, M'). The functional form of 
P(M, M') depends on which window function is used: It has a constant value for the case 
of the sharp- k space filter, while it var ies with the mass scale for the cases of more realistic 
top-hat spherical and Gaussian filters (IBond et al.lll99ll ). 
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Very recently, Maggiore &: Riotto ( 2010a b.cj) have generalized the standard excursion 
set formalism by incorporating both of the non-Markovian nature of the random walks and 
the stochastic nature of the absorbing barrier to account for the complexities associated 
with the actual formation and practical identification of the bound objects. Here, the non- 
Markovian random walking process corresponds to the case in which all the random steps 
in the trajectory are mutually correlated. Technically speaking, the non-Markovian random 
walking process depicts the gro wth of the density field smoothed with either the top-hat 
spherical or the Gaussian filter (IBond et al .1 1 1 9 9 it ) . Meanwhile, the growth of the density 
field smoothed by the sharp- space filter (i.e., the top-hat filter in the Fourier space) corre- 
sponds to a Markovian random walking process where only the adjacent walks are mutually 
correlated. 



Th is generalized excursion set mass function formalism proposed by lMaggiore fc Riotto 
(j2010aU blJd) was greatly successful not only in reproducing the N-body results but also in 
explaining physic ally the scale dependence of t he barrier heights in the ellipsoidal collapse 
models (see also. ICorasaniti fc Achitouvll2011aU bl). Unfortunately, however, this generalized 
formalism with the stochastic barrier is not directly applicable to the supercluster mass 
function for the following reasons. First of all, this generalized excursion set formalism 
resorts to the numerically deriv e d pro bability density distributions of 8 C from the N-body 
simulations of iRobertson et al.l (120091 ) to speci fy the stochastic scatte r of the absorbing 
barrier. The results of the numerical work of IRobertson et al.l (120091 ). however, is valid 
only for the cases of bound objects like clusters but not for the case of marginally bound 
systems like superclusters. 

Furthermore the stochastic absorbing barrier model is incapable of capturing the es- 
sential feat ure of the anisotropic for mation of the real superclusters. In this generalized 



formalism (IMaggiore fc Riotto 



2010bl ). it is assumed that the more massi ye an object is, th e 



more closely its collapse process follows the spherical dynamics (see also. lBernardeaulll994l ). 
In consequence, the stochastic barrier height 5 C has the smallest standard deviation around 
the spherical average value of 1.686 in the high-mass limit. However, the su perclusters ob 
serve d in the real Universe have quite non-spherical elongated shapes (e.g., lEinasto et al. 



20111 ). Moreover, recent high- resolution N-body simulations have revealed that the merging 
of the clusters into t he superclusters occur preferentially a long the cosmic filaments in an 
anisotropic way (e.g.. lWray et al.ll2006l ; iLee fc Evrardll2007l ). In other words, the marginally 
bound superclusters do not follow the spherical collapse dynamics unlike the bound massive 
clusters. 



In the following subsections, introducing a new concept of multi- dimensional clustering 
barrier to account for the anisotropic formation of the superclusters, we modify the excursion 
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set formalism to evaluate the supercluster mass function in a fully analytical manner. 



2.2. A Modified Excursion Set Formalism for the Superclusters 



The Zel'dovich approximation ( iZel'dovichl 1 1 9 70h . as the first order Lagrangian pertur- 
bation theory, relates the evolved mass density p at a given Eulerian region, x, to the initial 
shear eigenvalues, {Aj}f =1 (in a decreasing order), as 



p[x(q)] = pUl =1 [l - K(q)}- 



(5) 



where {Aj}? =1 are the three eigenvalues of the the initial tidal shear (i.e., linear deformation) 
tensor defined as the second derivative of the perturbation potential smoothed on a given 
mass scale at the corresponding Lagrangian region, q(x). Each of the initial shear eigenval- 
ues, {Aj}f =1 , is a random field, mutually dep endent on one another, satisfying the following 



joint probability density (jDoroshkevichl 1 19 TOT ) . 



p(Ai, A 2 , A 3 ;cr) 



3375 



8^' 



exp 



a 



iIj^j(Aj — A,-) 



(6) 



where a = a(M) in Equation 



This Zel'dovich model provides the simplest and the most intuitive explanation for the 
formation of the cosmic web in the Universe. The large scale structures that constitute the 
cosmic web such as voids, sheets and filaments are destined to form in the Lagrangian re gions 
where the three shear eigenvalues, {Aj}f =1 , have different signs (e.g.. lHahn et al.ll2007l ): the 
voids from the regions with all negative eigenvalues (Ai < 0), the sheets from the regions 
with two negative and one positive eigenvalues (Ai > 0, A2 < 0); the filaments from the 
regions with one negative and two positive eigenvalues (A2 > 0, A3 < 0). The gravitational 
collapse occurs along the directions of the eigenvectors (i.e., the principal axes of the linear 
deformation tensors) corresponding to the positive eigenvalues. The densest Lagrangian 
regions with three positive eigenvalues (A3 > 0) would condense out clumpy objects like 
clusters in the knobs of the cosmic web after collapsing along the three principal axes. 

The Zel'dovich approximation cannot be used to model the formation of the bound clus- 
ters due t o its notorious failure after th e first shell-crossing corresponding to the moment 
of Ai = 1 (IShandarin fc Zeldovichlll989l ). but it may be still adequately valid to depict the 
formation of the marginally bound superclusters. We assume that a supercluster of mass M 
forms in a Lagrangian region q where the following three conditions are satisfied: 
I. The three shear eigenvalues are all positive: A3 > A3 C = 0. 
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II. The largest shear eigenvalue does not exceed unity: Ai < Ai c = 1. 

III. The second to the largest shear eigenvalue exceeds one-half: A2 > \2c — 1/2. 

The first condition guarantees that a supercluster is not some section of a sheet or of a fila- 
ment but a three dimensional marginally bound object that are still under triaxial collapsing 
process. The second condition assures that the gravitational collapse along the major prin- 
cipal axis (the eigenvector corresponding to the second to the largest eigenvalue) has yet to 
be completed and thus the formation of a supercluster corresponds to an on-going merging 
event before the occurrence of the first shell-crossing. The third condition is derived not 
from the dynamical consideration of t he superclusters but from the consideration of their 



geometrical shapes. ISheth et al.l ( 120011 ) have shown that the prolateness, p, of a proto-halo 
region defined as p = (Ai — 2A2 + A3) / (25) is most likely to have zero value on average. Since 
the superclusters correspond to proto-halo regions, expected to end up as bound halos in 
the future, their prolateness should be most likely to be zero, too. Given that when the 
formation of a supercluster occurs when Ai c = 1 and A3 C = 0, the mostly likely zero value of 
the prolateness implies A2 C = 1/2. 

Now, let us translate the above three conditions for the supercluster formation into the 
excursion set jargons. The growth of the three shear eigenvalues is now described as multi- 
dimensional random walking process as a function of the monotonically increasing variable 
of a(M). A supercluster of mass M forms when the multi-dimensional random walk enters a 
clustering zone surrounded by two flat absorbing barriers (A2 C = 1/2 and A3 C = 0) and one flat 
reflecting barrier (Ai c = 1). Here a flat barrier means that its height is constant over er(M), 
while an absorbing (a reflecting) barrier means that it provides the lower (upper) limit on 
the corresponding shear eigenvalues. Figured] illustrates the entrance of a mult i- dimensional 
random- walk into the clustering zone. 

In our extended excursion set picture, the mass function of the superclusters, dN/dM, 
is directly proportional to the number density of those random-walks whose first entrance 
into the clustering zone occurs at a(M). The Jedamzik's excursion set formalism for the 
supercluster mass function is now written as 

f°° M' dN 

F(X lc , A 2c , A 3c ; M) = / dM'——P{M, M ) . (7) 

J M p dM 1 

Here -F(Ai c , A2 C , A3 C , M) is the cumulative probability of having the multi-dimensional random 
walk in the clustering zone at c(M), calculated as 

F(A lc ,A 2c ,A 3 ;M) = f " d\ x ! ' d\ 2 ! ' dA 3 p(A x , A 2 , A 3 ; a), (8) 

while P{M, M') is the conditional probability that a multidimensional random walk enters 
the clustering zone at some epoch a(M) after it touched the boundary of the clustering zone 
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at some earlier epoch of a(M') < cr(M), written as 



P(M,M') = P(A 1 >A lc ,A 2 >A 2c ,A 3 = 0|A , 1 = A lc ,A / 2 = A 2c ,A , 3 = 0) (9) 
" d\ 1 [ ^ d\ 2 [ ^ d\ 3 , p(Ai, A 2 , A3IA; = A lc , A' 2 = A 2c , A 3 = 0), (10) 

>X 2 JX 2c JO 

fx" dX i IxL C dX * Jo" rfA 3, p(Ai, A 2 , A 3 , A; = Aic, A' 2 = A 2c , A 3 = 0) 



p(Ai = Aie, A 2 = A 2c , A' 3 = 0) 



(11) 



where {Aj}f =1 and {A^}f =1 represent the shear eigenvalues on the mass scale of M and of M', 
respectively. 



Owing to the ingenious previous works of iDesjacquesI (120081 ) and iDesjacques h Smith 



( 120081 ). we have the analytic formula for the joint probability density distribution of the 
shear eigenvalues on two different mass scales, M and M' as, 



p(Ai, A 2 , A3, A x , A 2 , A 3 ) 



15 e 



(l- 7 2 )^(/3e_,ev,eA)e- 



3207T 2 (XV 6 



Here 



1 



7 : 



2vr W J 



xin^^-A.on^^-A'. 



d(ln k)k 3 P(k)W(k; M)W{k; M') 



Q 



157 



2(1 -I 2 
3 



4(1 -r») 



{5tr(A' 2 ) + 5tr(A' 2 ) - (trA') 2 - (trA) 2 + 2 7 (trA')(trA) } 



w(/3e_,6y,6 X ) 



2tt 



jo 



la 



dr I dtpexp 



30e. 



3/3e_e A 



(12) 

(13) 
(14) 
(15) 

(16) 



(with g = l + r 2 + e A '(l-r 2 )cos(2^), h = g 2 - 4(1 - e{,)r 2 ), 

where a' = a(M>), a = a(M), trA' = E;(A'/0, trA = EjA/a), tr(A' 2 ) = E* (A'/a') 2 , 
tr(A 2 ) = E,(A/^) 2 , e + = §(trA')(trA), e = §(trA' - 3A 3 /<x')(trA - 3A 3 /a), e v = (K - 
A 2 )/(trA' - 3A 3 /a'), ex = (Ai - A 2 )/(trA - 3A 3 /a). 

Figure [2] plots the conditional probability as a function of a(M) at 2 = for the five 
different cases of (Ai c , A 2c ). For all five cases, we let M' = 10 16 /i _1 M Q and A 3c = 0. As 
can be seen, the shape of P(M, M') changes sensitively with the mult i- dimensional barrier 
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heights. The random walks tend to enter the clustering zone at earlier epochs (i.e., smaller 
cr(M)) for the cases of the higher multi-dimensional barrier. It will be shown in §3] that the 
mult i- dimensional barrier heights are independent of redshift and background cosmology. 

For the practical calculation, we r eexpress Equatio ns © into a discrete form and regard 
it as a linear algebraic equation, as in iJedamzik! (119951 ): 



F(\ lc ,\ 2c ,\ 3c ;M i ) = P(M i ,M'] 



M' dN{M' 



dM> 



" dW- 



(17) 



where the Einstein summation rule is assumed over the subscript j. By inverting the matrix, 
P(Mj, M'-)M'-/S.M'-/p, we solve Equation (|T7|) to find the supercluster mass function, dN/dM. 



3. NUMERICAL TESTS 

To test the validity of our analytical model for the supercluster mass function, we com 
pare i t with the numerical results from th e three differ ent N-body simu lations: the Milieu 



nium flSpringel et all 120051 ). the CoDECS (1Baldill2011df ) and the MICE flCrocce et all 120101 ) 



simulations, all of which ran for the ACDM cosmology with slightly different cosmological 
parameters!!! TablefUlists the linear size of the simulation box, mass resolution, total number 
of particles, values of the three key cosmological parameters, and the halo-finding algorithm 
used for the three simulations. The Friends-of-Friends (FoF) algorithm with linking length 
of 0.2l p (where l p is the mean particle separation) was used for all of the three simulations 
to find the FoF groups of dark matter particles in which the bound halos are located. See 
the above three literatures for the more detailed descriptions of the simulations. 

The publicly available halo catalogs from each simulation provide such information on 
the resolved halos as their mass, positions, velocities and so forth at various redshifts. While 
the mass of each halo provided in the CoDECS and MICE catalogs corresponds to the FoF 
mass defined as the sum of the masses of all the particles belonging to the FoF group, the 
Millennium Run catalog contains not only the FoF mass but also the spherical overdensity 
(SO) mass of each halo defined as the mass enclosed by the spherical radius -R200 at which 
the halo mass density equals 200 times the mean background density at a given redshift. 

We conduct the following analyses using the catalogs from each simulation at three 
different redshifts (z = 0, 0.5, 1 for the cases of the Millennium and MICE catalogs while 



1 As for the CoDECS, the simulations ran not only for the ACDM cosmology but also for various coupled 
dark energy models. Here we use only the simulation data for the ACDM model. 
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z — 0, 0.44, 1 for the case of the CoDECS catalog). To begin with, we construct a mass- 
limited sample of the clusters from each catalog, selecting only those halos whose masses 
exceed the threshold value, M ct h = 1O 13 /i _1 M , including group-size and cluster-size halos 
but excluding the galactic halos for the supercluster membership. For the case of the MICE 
halo catalogs, however, a slightly higher threshold value, M Cjt h = 3.4 x 1O 13 /i _1 M , has to 
be used since all the halos of the MICE catalogs have masses larger than this value. 

We identify the superclusters as the clusters of clusters by applying the FoF algorithm 
to the clusters in each mass-limited sample at each redshift. The linking length of the 
FoF algorithm for the supercluster identification is set at l c /3 where l c is the mean cluster 
separation, following the conventional c r iterion suggested in the previous literatures (e.g., 



Kasun fc Evrardl 120051 ; IWray et al.ll2006t iLee fc Evrardl 120071 ) . The mass of each identified 



supercluster, M, is measured as the sum of the masses of its member clusters. The total 
number and mean mass of the superclusters from the three simulations at three different 
redshifts are listed in Tables EHU 

Binning In M and counting the number of the superclusters belonging to each differential 
mass bin, [InM, InM + cHnM], we determine the number density of the superclusters 
per unit volume, dN/dlnM, as a function of M at each redshift. To estimate the errors 
associated with the measurement of dN/dlnM, we also perform the Jackknife analysis: 
Dividing the superclusters into eight Jackknife subsamples, we determine dN/dhiM for 
each Jackknife subsample and calculate the one standard deviation scatter among the eight 
Jackknife subsamples as the errors associated with the measurement of dN/ d In M. 

Now, we would like to compare the numerical result from each simulation with the 
analytic model for the supercluster mass function. For this comparison, the analytic model is 
evaluated through Equations (IT|)- (TTT|) with the same cosmological parameters that were used 
for the given simulation and then renormalized to satisfy f?£ din M (dN/dlnM) = N tot /V 
where iV to t is the total number of the superclusters found in the numerical result and V is 
the simulation volume. 

Figure Elplots the numerical result of the supercluster mass function from the Millennium 
simulations (solid dots) and compare them with the analytic models (solid line) presented in 
§|2] at z = 0, 0.5 and 1 in the left, middle and right panel, respectively. The barrier heights in 
Equations (19|)-( ITT|) are set at Ai c = 1, \2 C = 0.5 and A3 C = at all three redshifts. Although 
the previous section (§|2]) focuses on the analytic model at z = 0, the supercluster mass 
function at higher redshifts can be readily evaluated just by replacing a(M) by a(M, z) = 
D(z)cr(M) where D(z) is the linear growth factor normalized to be unity at z = 0. The 



functional form of D(z) for the ACDM cosmology is given in lLahav et al.l ( I199ll ). As can 



be seen in Figure [31 the analytic mass functions of the superclusters agree excellently with 
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the numerical results at all of the three redshifts. Due to the relatively small box size of the 
Millennium simulations (Tabled]), however, the numerical results in the high-mass section 
(M > 10 15 h~ l M & ) suffer from large uncertainties, and thus it is hard to judge whether 
or not the analytic model matches the numerical results in the high-mass section for the 
Millennium case. 

FigureHlplots the same as Figure[3]but for the numerical results from the larger CoDECS 
simulations. The same barrier heights, A ic = 1, A 2c = 0.5 and A 3c = 0, are consistently 
implemented into the analytic model, while the key cosmological parameters are set at the 
values used for the CoDECS simulations (see Table[T]). The analytic models at three redshifts 
show excellent agreements with the numerical results from the CoDECS simulations, too. It 
is worth mentioning here that the value of the power spectrum amplitude is different between 
the two simulations as shown in Table [TJ o" 8 = 0.9 for the Millennium while <Tg = 0.809 for 
the CoDECS simulations. The fact that the analytic models with the identical thresholds, 
Aic = 1, A 2c = 0.5 and A3 C = agree very well with the numerical results from both of the 
simulations at three different redshifts indicates that the mult i- dimensional barrier heights 
determined from the Zel'dovich model is indeed independent of background cosmology and 
redshift. 

Figure [5] plots the same as Figure [3] but for the numerical results from the MICE 
simulations. Recall that a higher value of the mass threshold, M cX h = 3.4 x 10 13 h~ 1 M & , is 
used to construct the mass-limited sample of the clusters from the MICE simulations. This 
higher value of M c t h results in increasing the mean cluster separation, l c , which might in turn 
affect the multi-dimensional barrier heights. With the help of the ^-minimization method, 
the lower barrier heights of Ai c = 0.9 and A 2c = 0.45 are found to give the best agreement 
between the analytic model and the numerical results from the MICE simulations. Figure [5] 
reveals that our analytic models for the supercluster mass function agree impressively well 
with the numerical results from the MICE simulations, too, even in the high-mass section 
(M > 3 x 10 15 h^Mo) at all the three redshifts. 

The lower barrier heights for the MICE case can be understood by the following logic. 
The mass-limited cluster sample from the MICE simulation contains less number of the low- 
mass clusters (i.e., the group-size clusters) than those from the other two simulations and 
thus it has a larger value of the mean cluster separation, l c . In other words, the identified 
superclusters via the FoF algorithm with the fixed linking length of l c /3 must be less clus- 
tered. In the excursion set context, it can be said that their clustering zone is located in the 
lower area in the multi-dimensional space spanned by the three shear eigenvalues and cr(M). 
In consequence the random walks of the three eigenvalues are likely to enter the clustering 
zone in a relatively short <r(M)-interval. 
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We examine if the same higher value of M Cjt h is applied to construct the mass-limited 
cluster samples from the Millennium and CoDECS simulations, then the analytic model 
with the same lowered heights of the multi-dimensional barrier agree well with the numerical 
results. Figure E] plots the re-derived numerical results with higher cluster mass threshold 
of M C| th = 3.4 x 10 13 /i _1 M Q and compare them with the analytic model with lower barrier 
heights of Ai c = 0.9 and \2 C = 0.45. The value of A3 C is fixed at zero. As can be seen, 
the analytic model with lower barrier heights indeed agree very well with the numerical 
results from the Millennium and CoDECS simulations obtained using the higher cluster 
mass threshold. 

This result clearly shows that the mult i- dimensional barrier heights depend only on 
the strength of the clustering among the member clusters in the superclusters, revealing 
an intriguing aspect of our analytic model: the analytic model with the barrier heights of 
Aic = 1, A2c = 0.5, X^c = 0, derived purely from the Zel'dovich model match the numerical 
results obtained by including all the group-size low-mass clusters with mass threshold of 
10 13 h~ 1 M Q for the supercluster membership. 



4. RELATIVE ABUNDANCE OF THE RICH SUPERCLUSTERS 

As mentioned in given that the abundance of the rich superclusters in the Universe 
reflects how fast the structures grow and how frequently the clusters merge on the largest 
scale of the Universe, it may be possible in principle to use the relative abundance of the rich 
superclusters as a cosmological test of dark energy and gravity. To explore this possibility, 
we first define the relative abundance of the rich superclusters as 

r»r / w \ N(M>M sc ,z) 1 f°° „dN(M,z) , . 

5N lich (> M sc , z) = 1 ~ 5C; ' = — -j-r / dM y M \ 18 

AtotO) N tot (z) J Msc dM 

where 5A r i C h is the ratio of the cumulative mass function of the superclusters with masses 
larger than M sc to the total number of the superclusters, A to t- Note that the relative abun- 
dance of the rich superclusters is free from the renormalization of the supercluster mass 
function. 

Using the analytic model for the supercluster mass function derived in we can predict 
analytically <5A r i C h. We first investigate the variation of <5A r i C h with the key cosmological 
parameters in the fiducial ACDM model. Figure [7] plots the relative abundance of the rich 
superclusters SN rich at z = as a function of M sc (> 10 15 /i _1 M ) for four different cases of the 
density parameter Q m (left panel) and for four different cases of the linear power spectrum 
normalization factor as (right panel). For the evaluation of the analytic supercluster mass 



13 



funct ion, the other cosmological parameters are set at the WMAP7 values (IKomatsu et al. 
20 111 ) while the Zel'dovich condition of Ai c = 1, A 2c = 0.5, A 3c = are consistently used for 
the heights of the multi-dimensional barrier. As can be seen, the variations of o"g and fl m 
affect significantly the relative abundance of the rich superclusters SN^^ especially in the 
high mass section. As each of Q m and as increases, SN^ drops less rapidly with M sc . This 
result is consistent with the picture that the cosmic web grows faster and the clusters merge 
more frequently in a universe where the initial density fluctuations have higher amplitude 
and dark matter are more dominant. 

Assuming that our analytic model for the supercluster mass function also works in 
QCDM models (Quintessence+CDM) in which the Q uintessence scalar fie ld as a dynamical 
dark energy is responsible for the cosmic acceleration (jCaldwell et al.lll998l ). we also perform 
a feasibility study on how well SN ric ^ can constrain the dark energy equation of state. For this 
test, we focus on a particular toy QCDM model in which the dark energy equation of state 
is give n as w(z) = Wo + w-\z/(l + z) 2 where th e two parameters, wo and w\, have constant 
values (IChevallier fc Polarskill2001tlLinderll2003l ). For the evaluation of the supercluster mass 
function for this toy QCDM model, we use th e approxi mate analytic f ormula for the QCDM 
linear growth factor given in iBasilakosI ( 120031 ) (see also iPercivall 120051 ) . 



Figure |8] plots the relative abundance of the rich superclusters 5N Tic h at z = 0.5 for five 
different cases of the dark energy equation of states. The Jackknife errors from the MICE 
simulations are also overlapped with the analytic models to show explicitly how tight the 
constraints from the relative abundance of the rich superclusters would become if the same 
number of the superclusters were observed in the Universe. As can be seen, <5iV r i C h shows an 
appreciable change with w(z) in the high-mass section (M sc > 10 15 /i- 1 M Q ). As the value 
of Wi varies from —0.67 to 0.67, the value of SN Tich decreases by a factor of two at the mass 
scale of M sc = 3 x 10 15 h^M®. 

Finally, we also study how 5A" r ich changes in a toy modified gravity ( MG) model in 



which scales as a power l aw of the density parameter, D(z) oc fl^ (e.g., see iLinderl 12003 



20051 ; IShapiro et al.ll2010[ ). This toy MG model is distinguishable from the fiducial model 
(GR+ACDM) only in the value of 7 since the former has 7 = 0.68 while in the latter it 
is approximately 7 = 0.55. Assuming that our analytic model for the supercluster mass 
function also works in this toy MG model, we evaluate SN^^, which is plotted in Figure 
|9] (dashed line). The fiducial (GR+ACDM with WMAP7 parameters) case is also plotted 
with the Jackknife errors for comparison (solid lines) at z = 0.12 and z = 0.5 in the left 
and right panel, respectively. The Jackknife errors z = 0.12 and z = 0.5 are obtained from 
the CoDECS and MICE simulations, respectively. The difference in 5Arich between the two 
models at M sc = 3 x 10 15 /i _1 M reaches up to 50% and 66% at z — 0.12 and 0.5, respectively, 



-14- 



which proves clearly that the relative abundance of the rich superclusters at a given epoch 
is useful in principle as a cosmological test of gravity. 



DISCUSSION AND CONCLUSION 



We have developed a fully analytic model for the supercluster mass function by ex- 
tending the excursion set formalism to the case of marginally bound systems under triaxial 
collapsing process. The analytically derived supercluster mass function has been found to 
agree impressively well with three different N-body results at various redshifts, revealing its 
robustness. Using our model, we have shown that the relative abundance of the rich su- 
perclusters at a given epoch is sensitive to modification of gravity as well as to dark energy 
equation of state. 

We expect a wide application of the supercluster mass function in various fields. The 
abundance of the rich superclusters may be useful in constraining primordial non-Gaussianity 
parameter as we ll as neutrino mass due to its sensitivity to the initial cosmological conditions. 
As mentioned in lOguri et all (I2004T ). we also expect it to be useful in estimating the statistical 
significance of the amount of the warm hot intergalactic media in the superclusters which are 



the sources of soft X-ray and Sunyaev-ZePdovich effect (IMyers et al.ll2004t IZappacosta et al. 



20051 ; ISadeh &: Rephaeh1l2005l ). Our analytic model may also allow us to analytically evaluat e 
the late-time integrated Sachs- Wolfe (ISW) effect of superclusters (e.g.. lGranett et al.ll2008l ). 



Before applying our analytic model to the real Universe, however,we will have to under- 
take a couple of follow-up tasks. The first task is to investigate whether or not our model is 
valid for non-standard cosmologies. Given that our model is fully analytically derived from 
first principles, we have expected naturally it to work in QCDM or in models with modified 
gravity. However, this assumption has to be verified quantitatively using N-body tests. For 
instance, the Zel'dovich approximation in non-standard cosmologies might be different from 
Equation fl5J), which could lead to different conditions for the formation of the superclusters. 
The second task is to account for the projection effect along the directions of the line-of- 
sight on the mass measurement of the superclusters. Since most of the superclusters have 
elongated shapes along the cosmic filaments, it would be harder to find the member clusters 
due to the projection effect if a supercluster happens to be elongated along the direction of 
the line-of- sight. We plan to conduct these two follow-up works and to report the results 
elsewhere in the future. 



We acknowledge the use of data from the Millennium, CoDECS and MICE simulations 
that are publicly available at |http: / /www.millennium.com[ 
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http://www.marcobaldi.it/web/CoDECS.html, and http://www.ice.cat/mice, respectively. 



The Millennium Simulation analyzed in this paper was carried out by the Virgo Super- 
computing Consortium at the Computing Center of the Max-Planck Society in Garching, 
Germany. We also acknowledges the financial support from the National Research Foun- 
dation of Korea (NRF) grant funded by the Korea government (MEST, No.2011-0007819) 
and from the National Research Foundation of Korea to the Center for Galaxy Evolution 
Research. 
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Fig. 1. — Illustration of the mult i- dimensional random- walking process. The grey region 
corresponds to the clustering zone in the two dimensional projected space spanned by Ai 
and A2. 
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Fig. 2. — Conditional probability, P(M,M') in Equation ([7]), for the four different cases of 
the multi-dimensional barrier heights with A3 C = 0. 
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Fig. 3. — Comparison of the analytic mass function of superclusters (solid line) with the 
numerical results from the Millennium simulations (dots) at three different redshifts. The 
heights of the multi-dimensional barrier are given as Ai c = 1, A 2c = 0.5, A 3c = according 
to the Zel'dovich approximation. In each panel the errors represent the one standard de- 
viation scatter among eight Jackknife resamples. The mass threshold for the supercluster 
membership, M Cjt h is set at 10 13 /i _1 M Q . 



-23 - 




Fig. 4. — Same as Figure [3] but with the numerical results from the CoDECS simulations. 
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Fig. 5. — Same as Figure |3] but for the numerical results from the MICE simulations. The 
mass threshold, M C)t h, for the supercluster membership is 3.4 x 10 13 h~ 1 M Q from the MICE 
sample. The best-fit heights of the multi-dimensional barrier, Ai c = 0.9, A 2c = 0.45, A 3c = 0, 
are implemented into the analytic model of the supercluster mass function (solid line). 
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Fig. 6. — Same as Figures ISH but for the case that M Cjt h = 3.4 x 10 13 /i _1 M , the same value 
used for the MICE sample, is used for the construction of the mass-limited cluster samples 
from the Millennium simulation (left panel) and from the CoDECS simulations (right panel). 
The lower multi-dimensional barrier heights, Ai c = 0.9, \2 C = 0.45, A3 C = 0, that are used 
in Figure [5] are consistently implemented into the analytic models (solid line). 
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Fig. 7. — Relative abundance of the rich superclusters at z = for the four different cases 
of o"s (left panel) and for the four different cases of Q m (right panel). 
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Fig. 8. — Relative abundance of the rich superclusters at z — 0.5 for the four different cases 
of the dark energy equation of state, w(z) = w + W\zj{\ + z) 2 , assuming a toy QCDM 
model. The fiducial ACDM model (with wq = 0, w\ = and WMAP7 parameters) is also 
shown (solid line) for comparison. The errors are estimated as one standard deviation among 
8 Jackknife resamples from the MICE datasets. 
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Fig. 9. — Relative abundances of the rich superclusters for the case of the fiducial model 
(GR+ACDM with WMAP7 parameters) (solid line) and for the case of a toy model with 
modified gravity with D(z) = fi„ 68 (dashed line) at z — 0.12 and z = 0.5 in the left and 
right panel, respectively. The errors are estimated as one standard deviation scatter among 
8 Jackknife resamples from the CoDECS (left) and MICE (right) datasets. 
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Table 1. Simulation, linear box size, total number of dark matter particles, mass 
resolution, cosmological parameters and halo-finding algorithm. 



Simulation 


-^box 




M p 




n) 


Halo-Finder 




[/i _1 Mpc] 




[10 s h- l M & ] 








Millennium 


500 


10 10 


8.6 


(0.25, 0.9 


, 1) 


SO/FoF 


Codecs 


1000 


2 x 1024 3 


500.84 


(0.271, 0.809 : 


0.966) 


FoF 


MICE 


3072 


2048 3 


2300.42 


(0.25, 0.8, 


0.95) 


FoF 
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Table 2. Mass-limited sample, simulation, cluster mass threshold, total number of 
superclusters and mean supercluster mass at z — 0. 



sample 


simulation 


M C;th 
[W*h~ l M & ] 


iVtot 


Msc 

[l^h~ l M & ] 


I 


Millennium 


1.0 


35422 


5.7 


II 


Codecs 


1.0 


412337 


6.7 


III 


MICE 


3.4 


1851534 


13.7 


IV 


Millennium 


3.4 


9749 


14.0 


V 


Codecs 


3.4 


82016 


12.6 
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Table 3. same as Table 2 but for z = 0.5 (z = 0.44 for CoDECS). 



sample 


simulation 


M Cjth 

[io 13 ^-^©] 


N tot 


[lO 13 / i " 1 M ] 


I 


Millennium 


1.0 


26333 


4.5 


II 


Codecs 


1.0 


333739 


5.5 


III 


MICE 


3.4 


1091086 


11.3 
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Table 4. same as Table 2 but for z — 1. 



sample 


simulation 


M Cith 
[1O 13 /i~ 1 M ] 


iVtot 


[10 13 h^Mo] 


I 


Millennium 


1.0 


17518 


3.8 


II 


Codecs 


1.0 


214386 


4.4 


III 


MICE 


3.4 


494388 


9.3 



